home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / chebspec_r < prev    next >
Encoding:
Text File  |  1994-09-24  |  2.2 KB  |  80 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Chebyshev spectral differentiation matrix.
  4.  
  5. // Syntax:    C = chebspec ( N , K )
  6.  
  7. // Description:
  8.  
  9. //    C is a Chebyshev spectral differentiation matrix of order N.  
  10. //            K = 0 (the default) or 1.
  11. //          For K = 0 (`no boundary conditions'), C is nilpotent, with
  12. //              C^N = 0 and it has the null vector ONES(N,1).
  13. //              C is similar to a Jordan block of size N with eigenvalue zero.
  14. //          For K = 1, C is nonsingular and well-conditioned, and its eigenvalues
  15. //              have negative real parts.
  16. //          For both K, the computed eigenvector matrix X from EIG is
  17. //              ill-conditioned (MESH(REAL(X)) is interesting).
  18.  
  19. //      References:
  20. //      C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral
  21. //             Methods in Fluid Dynamics, Springer-Verlag, Berlin, 1988; p. 69.
  22. //      L.N. Trefethen and M.R. Trummer, An instability phenomenon in
  23. //             spectral methods, SIAM J. Numer. Anal., 24 (1987), pp. 1008-1023.
  24. //      D. Funaro, Computing the inverse of the Chebyshev collocation
  25. //             derivative, SIAM J. Sci. Stat. Comput., 9 (1988), pp. 1050-1057.
  26.  
  27. //    This file is a translation of chebspec.m from version 2.0 of
  28. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  29. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  30.  
  31. //-------------------------------------------------------------------//
  32.  
  33. chebspec = function ( n , k )
  34. {
  35.   local (n, k)
  36.   global (pi)
  37.  
  38.   if (!exist (k)) { k = 0; }
  39.  
  40.   // k = 1 case obtained from k = 0 case with one bigger n.
  41.   if (k == 1) { n = n + 1; }
  42.  
  43.   n = n-1;
  44.   C = zeros(n+1,n+1);
  45.  
  46.   one = ones(n+1,1);
  47.   x = cos( (0:n)' * (pi/n) );
  48.   d = ones(n+1,1); 
  49.   d[1] = 2;
  50.   d[n+1] = 2;
  51.  
  52.   // eye(size(C)) on next line avoids div by zero.
  53.   C = (d * (one./d)') ./ (x*one'-one*x' + eye(size(C)));
  54.  
  55.   //  Now fix diagonal and signs.
  56.  
  57.   C[1;1] = (2*n^2+1)/6;
  58.   for (i in 2:n+1)
  59.   {
  60.     if (mod(i,2) == 0)
  61.     {
  62.       C[;i] = -C[;i];
  63.       C[i;] = -C[i;];
  64.     }
  65.     if (i < n+1)
  66.     {
  67.        C[i;i] = -x[i]/(2*(1-x[i]^2));
  68.     else
  69.        C[n+1;n+1] = -C[1;1];
  70.     }
  71.   }
  72.  
  73.   if (k == 1)
  74.   {
  75.     C = C[2:n+1;2:n+1];
  76.   }
  77.  
  78.   return C;
  79. };
  80.